%% main_householdbalancesheet.m
%
% Laibson, Maxted, and Moll
% Main code: "Present Bias Amplifies the Household Balance-Sheet Channels of Macroeconomic Policy"
%--------------------------------------------------------------------------

clear;

%--------------------------------------------------------------------------
%% Model Set-Up

%%%%%%%%%%%%%%%%%%%%%%
saveOutput = 1;

calibrationCase = 1;  %see model_calib.m for details
sim_case = 0.0;       %see initialize_sim_case.m for details
%%%%%%%%%%%%%%%%%%%%%%

addpath('subfunctions')
model_calib
initialize_sim_case
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
%% Solve HJB Equation and Steady State Distribution
if HJB_preSolved ~= 1 %Skip this block of code if already solved

    solveHJB;
    interventionMatrix;

    %----------------------------------------------------------------------
    % KF (steady state)
    rb_spot_ind = rb_spot_ind_stationary;
        %Create large transition matrix A including rm switches, given r-spot
        A = sparse([]);
        for rm_ind = 1:Nr
            A = [A; [sparse(M,M*(rm_ind-1)), AStore{rb_spot_ind,rm_ind}, sparse(M,M*(Nr-rm_ind))]];
        end
    
    deathFlag = 1; %Households retire at rate la_BY
    if useAlternateStationary ~= 1
        findStationary = 1;
        solveKF
        g_stationary = g;
    else
        g_stationary = gAlternate;
        g = g_stationary;
    end

    %Calculate average net worth
    netWorth = sum(sum(sum(sum(((h - aaa)+bbb).*g*db*da))));
    disp(['Average Net Worth = ' num2str(netWorth)]);
    %----------------------------------------------------------------------
    
else
    %Keep current set-up
    saveOutputCurrent = saveOutput;
    sim_caseCurrent = sim_case;
    
    %Load in pre-solved model
    if ~exist('altLoadName', 'var')
        load(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll']);
    else
        load(['output/', altLoadName]); 
    end
    
    %Reset set-up details
    saveOutput = saveOutputCurrent; clear saveOutputCurrent;
    sim_case = sim_caseCurrent; clear sim_caseCurrent;
    initialize_sim_case
    
    if useAlternateStationary == 1
        g_stationary = gAlternate;
        g = g_stationary;
    end
end    

%Save steady-state workspace
AStoreSS = AStore;
interventionStoreSS = interventionStore;
cStoreSS = cStore;
sssStoreSS = sssStore;


%--------------------------------------------------------------------------
%% Simulate Economy Forward
simulateForward;


%--------------------------------------------------------------------------
%% Save Workspace
if saveOutput == 1
    close all
    if modelSaveSmall == 1
        saveList = {'meanc_t_store','meana_t_store','meanb_t_store','meannw_t_store','rm_t_store', ...
            'Delta_path','rb_spot_ind_path','time_path', 'g_stationary', 'g_t_store'};
        save(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll', sn_remainder_sim], saveList{:})
    else
        save(['output/beta', num2str(1000*beta), 'rho', num2str(10000*rho), '_outputAll', sn_remainder_sim])
    end
end    


